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Computational models at different space-time scales allow us to understand the 
fundamental mechanisms that govern neural processes and relate uniquely these 
processes to neuroscience data. In this work, we propose a novel neurocomputational 
unit (a mesoscopic model which tell us about the interaction between local cortical nodes 
in a large scale neural mass model) of bursters that qualitatively captures the complex 
dynamics exhibited by a full network of parabolic bursting neurons. We observe that the 
temporal dynamics and fluctuation of mean synaptic action term exhibits a high degree of 
correlation with the spike/burst activity of our population. With heterogeneity in the applied 
drive and mean synaptic coupling derived from fast excitatory synapse approximations 
we observe long term behavior in our population dynamics such as partial oscillations, 
incoherence, and synchrony. In order to understand the origin of multistability at the 
population level as a function of mean synaptic coupling and heterogeneity in the firing rate 
threshold we employ a simple generative model for parabolic bursting recently proposed 
by Ghosh et al. (2009). Further, we use here a mean coupling formulated for fast spiking 
neurons for our analysis of generic model. Stability analysis of this mean field network 
allow us to identify all the relevant network states found in the detailed biophysical model. 
We derive here analytically several boundary solutions, a result which holds for any number 
of spikes per burst. These findings illustrate the role of oscillations occurring at slow time 
scales (bursts) on the global behavior of the network. 
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1. INTRODUCTION 

The neuronal spike-burst activity is characterized by recurrent 
transitions between rest state and firing state where bursts are 
temporal groupings of multiple spikes. Certain cells in the mam- 
mal brain, for example, neurons in the thalamus during periods of 
drowsiness, attentiveness, and sleep are known to exhibit this type 
of spike-burst behavior (Sherman and Koch, 1986; Steriade and 
Llinas, 1988; McCormick and Feeser, 1990; Steriade et al, 1993; 
Amzica and Steriade, 1998). Autonomously bursting neurons 
are found in a variety of neural systems, from the mammalian 
cortex (Morris and Lecar, 1981; Dhamala et al., 2004a,b) to brain- 
stem (Hindmarsh and Rose, 1984; Wang, 1994; Izhikevich, 2007; 
Jirsa and Mcintosh, 2007; Jirsa, 2008). When neurons are cou- 
pled with each other, they produce different modes of behavior, 
including synchrony and phase-locking, which have been impli- 
cated in memory, cognition, sensory processing, motor planning, 
and execution (McCormick and Feeser, 1990; Wang, 1994; Jirsa 
and Mcintosh, 2007). Many neurological diseases, on the other 
hand, including Parkinson, schizophrenia, and epilepsy, are the 
result of abnormal synchronization (Uhlhaas and Singer, 2006; 
Jensen et al, 2007), which suggests that a better understand- 
ing of the basic mechanisms producing synchrony and phase 
locking will be a stepping stone toward the repair of brain func- 
tion. Modeling attempts using large scale networks to understand 



emergence of cognitive states rely heavily on the approximation 
of the dynamics as a neural ensemble. The concept of a neu- 
ral mass like abstraction (Hebb, 1949; Beurle, 1956) designates a 
group of Co-activated neurons capable of acting like a closed sys- 
tem when performing a certain function. A small scale network 
of this kind is sometimes referred to as a "neurocomputational 
unit." In large scale brain networks, these mesoscopic units of 
operation serve as the network nodes (see for instance, Deco 
et al., 2008, 2011; Ghosh et al., 2008). On intermediate spatial 
scales of few cm, neural activations along the spatially continuous 
cortical sheets are described by neural fields, for which the con- 
nectivity is assumed to be translationally invariant (see, Wilson 
and Cowan, 1972; Nunez, 1974; Amari, 1977; Jirsa and Haken, 
1997; Feng et al, 2006; Jirsa, 2009; Robinson, 2011). To define 
such small neurocomputational units, simplified neuron models, 
known as phase models, offer an attractive tool for the study of 
network modes, since they allow for detailed mathematical analy- 
sis of network dynamics (Breakspear et al., 2010). As an example, 
Carbal et al. have explored the role of local network oscillations 
in resting-state functional connectivity by using such phase oscil- 
lators in the respective nodes of the simulated network. They 
have shown when these oscillatory units are integrated in the net- 
work, they behave as weakly coupled oscillators. Moreover, for a 
set of network parameters they found subsets of nodes tend to 
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synchronize although the network is not globally synchronized 
(Cabral et al., 2011). For the present work we use a recently 
proposed phenomenological model that admits parabolic burst- 
ing in one dimension, which is a type of bursting observed 
in the R-15 neuron in abdominal ganglion of aquatic mollusc 
Aplysia Californica (Ermentrout and Kopell, 1986; Izhikevich, 
2000; Ghosh et al., 2009). This type of bursting can arise even 
without bistability in the generation of spikes. The investigation 
carried out in this work with a detailed neuron model capable of 
displaying spiking and bursting behavior and a minimal model 
that not only reproduces the mean field amplitude of the original 
networks but also capture the most important temporal features 
of its dynamics. The detailed model used here is extensively dis- 
cussed in Rinzel and Ermentrout (1989). On the other hand, our 
phase model is a minimal model that captures the generality of 
the mechanism of bursting present in the detailed model. As we 
vary network parameters including mean field coupling strength 
and dispersion, both networks display various temporal dynam- 
ics. In order to understand these states in mathematically tractable 
terms we take advantage of the mean field coupled network of 
phase model. Our goal is to identify to what degree this mean field 
model serves as a reliable neurocomputational unit and captures 
the qualitative features of temporal dynamics of the full network 
as a function of the investigated network parameters. Mean field 
analysis for singleton burst reveals solutions such as incoherence 
and partial oscillation which can be completely described ana- 
lytically. However, as we are interested in a multispike system 
where analytical calculation is rather non-trivial and therefore, 
we combine semi-analytical approach with numerics to derive 
the stability diagram. Mean field phase network allow us to iden- 
tify the mechanism of transitions between various network states 
that appear as solutions of the full network. Stability diagram 
is independent of number of spikes per burst and qualitatively 
commensurates well with the findings in our full network. The 
paper is structured as follows. In the next section, we introduce 
the Rinzel-Ermentrout model (Rinzel and Ermentrout, 1989) for 
parabolic bursting and describe the model in details. In the fol- 
lowing section, we couple individual neurons via global coupling 
and present our analysis of this network model. In the subse- 
quent section, we set up a generic network of bursters coupled to 
their mean field and derive semianalytically all the network states 
and corresponding phase transition boundaries. In the next sec- 
tion, we derive numerically a stability diagram using global phase 
coherence measure. In the final section, we summarize the results 
obtained from mean field descriptions and link them systemat- 
ically with the network states obtained from biophysical model 
network. 

2. MATERIALS AND METHODS 
2.1. SINGLE NEURON MODEL 

A dynamical system with multiple time scales (for example, a 
neuron with spiking-bursting behaviors) can be written in a sin- 
gularly perturbed form: x = f(x, y), y = rg(x, y), where x is the 
vector of fast variables, y the vector of slow variables that modu- 
late the fast activity, and r 1 is a ratio of fast/slow time scales. 
A system which has been proposed to describe parabolic burst- 
ing behavior is known as Rinzel model (1989). Single neuron 



model parameters used here are exactly as described in Rinzel and 
Ermentrout (1989). 

V = (I - lea ~ igKW + gkcaZKV - V K ) - gl (V - Vj))/c 
w = (|>Ooo - w)/x w 

Ca = €(-(jL7 Cfl - Ca) 
h = €(noo(V) - n)/t„ (1) 

where I Ca = (gCaMcotV) + &c««)(V - V Ca ), z = Ca C + Caa and 
gating functions are 

moo(V) = 0.5(1 +tanh((V-vi)/r 2 )) 
Woo(V) = 0.5(1 + tanh((V - v 3 )/v 4 )) 
Woo (V) = 0.5(1 + tan h((V- vs)/v 6 )) 
t w (V) = 0.5(1 + tanh((V - v 3 )/2v 2 )) (2) 

where V is the membrane potential, w is associated with the 
fast current, Na + or K + , Ca and n are the two slow currents, 
Model parameters which are held fixed throughout our sim- 
ulations are, V K = -84, V/ = 60, V Ca = 120, g K = 8, gi = 2, 
c = 20, Vi = 1.2, v 2 = 18, v 6 = 24, v 5 = 12, v 3 = 12, v 4 = 17.4, 
x„ = 0.05, 4> = 0.06666666, gCa = 4.0, |x = 0.025, Ca 0 = 1, € = 
0.0005, smdg kCa = l,g sCa = 1. 

/ is the applied input current. The ionic currents are given 
by an ohmic leak current, determined by the leak conductance 
gl and leak reversal potential V;, and a Na + current which is 
responsible for the generation of spikes. The dynamics of this 
model which is relevant to our study is outlined as follows. When 
the input current I exceeds a critical value I c a single neuron 
described by Equation (1) undergoes a Saddle-node bifurcation 
on an invariant circle (SNIC). This same system for two differ- 
ent parameterization of I and in the presence of the slow currents 
can exhibit both spiking as well as parabolic bursting behavior. 
Spiking behaviors are elicited for a slightly higher value of the 
external drive. For example, to observe a typical burst-like pat- 
tern in this system we held the input current to the values I = 68 
and for spikes I > 70. Figure 1 displays the relationship between 
the applied input current and a parabolic bursting pattern that is 
observed in the single neuron dynamics. 

2.2. PHASE MODEL 

The generality of the underlying mechanism for parabolic burst- 
ing is investigated in details by numerous authors (Ermentrout 
and Kopell, 1986; Baer et al, 1995; Izhikevich, 2000). In many 
such formulations, parabolic bursting neurons are typically in 
their canonical form described as: 

9 = [l-cos(9) +f(x,y)] 

x = ix x [x ); (e) -x] 

y= \y-y\y n ($) -y] O) 

where function /(x, y) in the above equation couples to spike gen- 
erative mechanism depending on the slow variables x, y dynamics, 
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respectively. The function / (x, y) is a smoothly varying peri- 
odic function alternating signs such that the system undergoes 
a SNIC to generate parabolic burst at the single neuron level. 
Recently Ghosh et al. (2009) has also proposed a simpler model 
that in principle captures the underlying mechanism of parabolic 
bursting involving only a circular phase variable G and more- 
over, involve only one slow term to allow the fast dynamics to 
enter or get out of repetitive firing. Motivation for using such a 
model is primarily mathematical tractability. Parameter space of 
this model cannot be directly linked to the biophysical parame- 
ters, however, qualitatively it may account for the transient and 
longterm behavior of more detailed biophysical models. In this 
model a single neuron is described by the following equation, 



7 — cos 9 ■ 



(4) 



In Equation (4) a slow variable activation term is represented by a 
modulation term cos(-) which mimics the entire slow subsystem 
instead of describing it as a separate dynamical system, 7 is the 
applied input current and n is an integer, which determines the 
number of spikes per burst. In our simulation with this model 
all the results are for n = 5 spikes per burst unless otherwise 
specified. 

2.3. FULL NETWORK MODEL 

Golomb and Rinzel (1993) considered a heterogeneous net- 
work of all-to-all coupled inhibitory bursting neurons and 
found regimes of synchronous, anti-synchronous and asyn- 
chronous behavior when the width of the heterogeneity was 
changed (Golomb and Rinzel, 1993; Stefanescu and Jirsa, 2008, 
2011; Smeal et al., 2010; Jirsa and Stefanescu, 2011). We describe 
our network equations via a fast instantaneous coupling. N 
synaptically coupled (all-to-all) parabolic bursting neurons are 
described by a similar set of non-linear differential equations with 
fast chemical synapse. To this end we formally describe: 

Vi = (bU - I Ca ~ (gKW + gkcaZ)(V, - V k ) 

-gliVi-Vd+KSWi-V&V/c 

Wi = tyiWoo - Wi)/X w 

Cat = €(-mJ m - Cfl,) 
n, = € (Moo(V) - «,) /x„ 

s, = a s (V,)(l-5i)-| (5) 
P 

where all the parameters and the gating variables inherit from 
the single neuron model Equation (1,2) and b is a rescaling fac- 
tor to applied drive to cross the threshold and elicit spike/burst 
in the uncoupled system. Stimulus that all the neurons see 7; > 
0 are drawn from a uniform distribution assumed to be sym- 
metrically distributed over the interval 7, e [2.1 — A7, 2.1 + Al]. 
Where A7 is the spread of the applied stimulus parameter. A7 
introduces a heterogeneity in the spike threshold. The synaptic 
coupling appears as an ensemble average given by S = ^ Ylf= l s »> 



The synaptic strength K is the same for all the neurons. For the 
entire simulation, we fixed the reversal potential of potassium 
ions to v t h ~ 0.0 (for purely excitatory connectivity). 

Analysis is carried out for a fast synapse (AMPA-type glu- 
tamate receptors), such as those found in the auditory system, 
the rise time is instantaneous, and post-synaptic responses com- 
mence almost instantaneously after the start of presynaptic action 
potential (Nunez, 1974; Morris and Lecar, 1981). This brisk com- 
munication is a consequence of rapid calcium-channel kinetics, 
which allows significant calcium entry during the upstroke of 
the presynaptic action potential (Sabatini and Regehr, 1996). 
Under the fast synapse approximation the variable s; relaxes 
much more rapidly than Vj, in which case we may apply a 
quasi-static approximation to (Equation 5) (e), s,- ~ 0, allow- 
ing us to adiabatically eliminate the synaptic variable via s, = 

(i + p + exp(-y-/2)) • ^ e ti me course °f tne postsynaptic conduc- 
tivity caused by an activation of AMPA receptors can be captured 
by a rise time p r ; se = 0.09 ms and decay time p 1 decay = 1.5 ms 
(Gabbiani et al., 1994; Parnas and Parnas, 1994). Numerical 
results in Figure 3 provides a good approximation for f$ in the 
range between [0.01ms, 0.5 ms]. Although, we have provided 
here the details about the fast excitatory synaptic connectivity, 
our approach can be readily extended to inhibitory connec- 
tivity as well. In the continuum limit, a mean field formula- 
tion with inhibitory synaptic coupling is provided in details in 
Appendix. 

2.4. MEAN FIELD COUPLED PHASE MODEL 

Each generic neuron is coupled to this mean field and typi- 
cally their response to the mean field expressed as R(Q) explicitly 
dependent on 6, and implicitly on time. In absence of any cou- 
pling, their vector field flow on a real line is governed by F(9) = 
w — cos(9) — cos(9)/n. In the absence of the term cos(9)/« 
expression reduces to a mathematical description used in Roy 
et al. (2011). Together, we can write for N (still finite) such 
neurons: 



9,=F(9,)-r7?(9,), 



(6) 



Recently, we have proposed a formulation for mean synaptic acti- 
vation term under fairly general setting and taking advantage of 
instantenous activation, deactivation between pre and postsynap- 
tic events. It allows one to describe synaptic activation variable 
5; = P ,_ v . as a non-linear transfer function of membrane 

l+P+exp(Vj 

voltage (Roy et al., 2011). Moreover, we have described how the 
mean field coupled spiking neurons can be described mathe- 
matically with this synaptic coupling. Details of this formulation 
is described elsewhere, (Roy et al., 2011). Collective activity of 
synapses is described by a mean field. For a given population of 
neurons is expressed more formally as, 



K N 



p 



N 



1=1 



(l + P + exp(-^)y 



i # I (7) 



where a $i (Vi) 



(l + exp(-Vj/2)) 



is a sigmoidal activation function. 



where T is the mean field influence function. Coupling K is 
the same for all the neurons. In our previous work, response to 
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such mean field coupling explicitly described as R(Q,) = sin G; 
(cos 6, - v th ), 

9, = F(Qi) - T sin 9,(cos G, - v th ) + O(e), (8) 

where 0(e) contains non-circular deviations of the order 6 that 
results due to perturbations. v tn ~ 0.0 for all simulations and 
analytical calculations unless mentioned otherwise. It is impor- 
tant to note that the couplings in the phase descriptions retain 
their mathematical expression in the full model plus some linearly 
added correction terms, which scale with the degree of order of 
deviation from the circle (Roy et al., 201 1). Hence, in application 
it is rather suitable when phase perturbations are close to the cir- 
cular orbit. The above equation further can be written combining 
the terms containing a single Fourier harmonic in the coupling 
plus the higher order Fourier terms. 

9,- = CO; — sin 9; — sin(9,/ «) + P(9/) sin 9, v tn 

+ 0(29,-) + 0(c), (9) 

K N ft 
w ,tt(l+P + exp(-^)) 

See for details (Roy et al., 2011). Where, in Equation (6) the 
frequencies T,- > 0 are assumed to be symmetrically distributed 
over the interval 7, e [7 — A7, 7 + AT] according to a uniform 
probability distributions. 

2.5. CHARACTERIZATION OF SPIKE/BURST COHERENCE IN 

BIOPHYSICAL NETWORK MODEL WITH MEAN FIELD COUPLING 

The bursting coherence and incoherence is quantitatively charac- 
terized in terms of a statistical-mechanical spike-based measure. 
We consider an excitatory population of neurons coupled to a 
common mean field drive and heterogeneity in their thresh- 
old for spikes/bursts. By varying the strength of the coupling 
K and the stimulus spread A7 we investigate the emergence of 
spike/burst coherence. Emergence of collective spiking/bursting 
coherence may be well described by the (population-averaged) 
global potential, 

1 N 

^mean(f) = - V V, (11) 

TV ^ 

>= 1 

In the thermodynamic limit (TV — > oo), a collective state becomes 
coherent if SV me3B (t) = [V mean (f) - V mean (f)] is non-stationary 
(i.e., an oscillating global potential Vmean appears for a coher- 
ent case), where the overbar represents the time average, and 
also, the correlated mean field r(r) activity appears oscillatory. 
Otherwise (i.e., when Vinean ls time independent or stationary), 
it either becomes incoherent (IN) or partial oscillatory (PO). In 
TV — > oo limit both these states converges to a stationary solu- 
tion. Thus, the mean square deviation of the global potential is 
a global marker for mean burst coherence for the entire pop- 
ulation described here. More formally one can write it as (i.e., 
time-averaged fluctuations of V m ean)> 

R(t) = (Vmean(f) - W*)) 2 (12) 



plays the role of an order parameter used for describing the 
coherence-incoherence transition (Manrubia et al., 2004). For the 
coherent (IN) state, the order parameter _R(r) approaches a non- 
zero (zero) limit value as N goes to the infinity. We compute R(t) 
in Equation (12) as a function of mean field coupling strength K 
and dispersion parameter AT for the full system. We vary both 
K, AT from 0 to 1 in a step size of 0.01. Subsequently, computed 
values of R(t) is plotted in grid size of 100 x 100. Contour plot 
is colorcoded from low values at zero (blue) to high values at 1 
(red). Nearly (in phase or anti phase) synchronized population 
spike/burst activity is lumped into a regime with labeled as SR and 
IN population spike/burst activity is lumped into a regime called 
IN activity. In the IN regime as described above R(t) values stays 
close to zero with substantial subthreshold fluctuations. Partial 
bursty regime is labeled as PO observed for R(t) values stationary 
and close to values other than zero. This regime displays dynami- 
cal behaviors far from synchrony, such as multi-clustering (some 
of the neurons are firing incoherently while others are not firing 
at all) in the phase for instance. Depending on the heterogeneity 
in stimulus spread we get random distribution of phases such that 
individual members can exhibit cluster hopping. Multiclustering 
in our model can reliably be captured using an ensemble average 
quantity rotation number p,- given by Equation (14). 

2.6. CHARACTERIZATION OF SPIKE/BURST COHERENCE IN PHASE 
NETWORK MODEL WITH MEAN FIELD COUPLING 

The bursting coherence and incoherence is quantitatively charac- 
terized in terms of statistical mechanical order parameter coher- 
ence measure. As an alternative to storing and plotting many 
time series data 9,(f), i = 1, . . . , TV for all TV = 1000 variables, we 
define an order parameter 

1 N 

TV t— 1 

>= l 

Equation (13) measures the population dynamics. The advan- 
tage of using such a formulation becomes apparent immedi- 
ately. Let's say our model system has periodic orbit then 9,(f) 
9,(f + T), where T periodic pacing spikes or bursts (latency). 
Then in order parameter space one can can detect this state in 
a straight forward manner as a solution 7?e(f)T?e(f + T). This 
result holds for all i, t. In this case, £9 dynamics is dominated 
mostly by the x co-ordinate dynamics. Absolute values of mean 
order parameter mod Rq < I. There is a mathematical rela- 
tionship of macroscopic global phase measure with macroscopic 
Vmean (0 in Equation (11). The interval between each micro- 
scopic spike/burst in an arbitrary £th stripe of spike/burst can be 
determined in a statistical-mechanical way by taking into con- 
sideration its contribution to the macroscopic global membrane 
potential V me an(£)- In this interpretation, the time series of the 
global potential Vmean (0 has a local maxima and minima, respec- 
tively and strictly bounded between [0,1]. The global cycle in the 
suprathreshold regime starting from the minimum of Vmean (£) 
which appears first after the transient time is regarded as the first 
global cycle, which is denoted by G\. The 2nd global cycle G2 
begins from the next following right minimum of Gl, and so on. 
Then, we can introduce an instantaneous global phase measure 
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6(f) of V mean (f) via a linear interpolation in the two successive 
subregions forming a complete global cycle (Lim and Kim, 201 1). 
A microscopic spike makes the most constructive (inphase) con- 
tribution to Vmean when the corresponding global phase Qk for 
fcth cycle of spikes/burst is 2nu (n = 0, 1, 2, . . .), while it makes 
the most destructive (anti-phase) contribution to V me3a (t) when 
6, for an arbitrary ith cycle of burst is 2(nl/2)jt. By averaging the 
contributions of all microscopic spikes within a burst in the ith 
burst stripe to Vmean > we can obtain the following degree of order- 
ing of spikes/bursts. Hence, the contribution of fcth microscopic 
burst occurring at the time tk is ordered by R^fa). If the degree of 
synchrony is high between the bursts/spikes then -Re(fjt) — > 1. We 
quantify the average firing frequency to compare the long-term 
behavior of individual neurons in the population model. We com- 
pute the average frequency (also known as the rotation number) 
of population of neurons using 

6; 

Pi = lim — , i = 1, . . . ,N. (14) 

t^oo t 

Averaging is carried out over about 1000 neurons starting from 
random initial conditions after the transient have died out. 
Collective states of ensemble of N = 1 000 neurons with spikes per 
burst n = 5 as indicated by their rotation numbers with uniform 
distribution of frequency I in the interval [2.1 — AI, 2.1 + AI]. 
Different branch of rotation index indicate different dynamical 
states of the network as a function of mean field coupling strength 
K, AI. We carry out a grid search in the 2D parameter space K, 
A7. Our goal is to obtain a phase transition diagram to under- 
stand long-term collective behavior of Equation (8) for large N, as 
a function of the coupling strength K > 0 and the stimulus spread 
AI € [0, 1). Global order parameter -Re(f) is computed for differ- 
ent parameterization of K, AI and embedded on a contour plot. 
Color spectrum is the same as the one used for displaying phase 
diagram in the full network. The values which are high and close 
to 1 are indicated by red and the values which are close to zero are 
indicated by blue. 

2.7. CLUSTERING ANALYSIS IN N COUPLED FULL AND PHASE 
NETWORK MODEL 

We describe firing patterns in large networks (finite N) with exci- 
tatory mean field coupling in terms of array diagrams. Array 
diagrams are obtained by simulating a coupled system consisting 
of mean field coupled biophysical neurons (N = 100) governed 
by the Equation (5). All the coupling coefficients are the same 
K where i = 1 , . . . , N. In the arrays the intensity of the voltage 
variables V\ , . . . , V; have been encoded in color spectrum. Two 
different color spectrums are used for the biophysical network 
(see Figure 4). In the first color spectrum blue part of the array 
values implies the quiescent activity of the spikes where the volt- 
age variables have relatively lower values. All the other colors in 
the spectrum indicates the higher values for the voltage variables, 
consequently these pixels in the array imply the spike activity. 
The horizontal line of the array shows the time with increasing 
epoches of activity. The second color spectrum used here shows 
burst depiction in the nearly coherent parameter regime. Green 
colors in the array indicate completely silent neurons. Purple pix- 
els on the green background shows burst activity. On the vertical 



axis neuron index are aligned and again, on the horizontal axis 
gives the direction of time. These diagrams were obtained from 
a phase network by monitoring phases of individual neurons 

i = 1 N and aligning them on the vertical axis. The choice 

of the color spectrum used for phases is given by a colorbar with 
uniformly distributed phase values. In Figure 9 red color index 
in the spectrum corresponds to higher phase values of 6 (close 
to jt) and orange color index are for lower phase values (close 
to — jt). First initial conditions 6,(0) is generated randomly and 
then they are sorted according to their neuron index and subse- 
quently distributed uniformly about [— jt, jt]. The parameters K, 
AI, for both realizations are chosen from SR, IN regime of the 
respective phase diagrams. 

2.8. NUMERICAL PROCEDURES AND VISUALIZATION OF THE SYSTEM 
DYNAMICS 

Two network models were implemented in Matlab, numerically 
integrated using second order Runge Kutta routine and Euler- 
Maruyama (EM) method (Higham, 2001). The simulations were 
performed with a fixed time step of dt = 0.05. The first 200 time 
points of the simulation are disregarded to set the network to 
a steady state. Thus, the results within this time were ignored. 
The membrane potential V(t), standard deviation of membrane 
potential std V, mean field T(f), order parameters R(t), R%(t) are 
captured for the entire population. For full network, simulation 
is carried out for N = 100 neurons and for the phase network 
for N = 1000 neurons. Numerical Phase diagrams are obtained 
using parallel for loops implemented in Matlab. Coupled mean 
field Phase model represented in Equation (8) can be visualized 
as a collection of N points rotating around the unit circle, where 
the estimated phase for each neuron 6,(f) denotes their position 
on a ring or a circle at time t. This alternate representation of 
the dynamical system (as N points moving along a circular refer- 
ence frame, instead of a single point tracing out a trajectory in an 
N-dimensional phase space) is possible because the system's state 
space, the N-torus, is equivalent to N copies of the unit circle. It 
is worth noting that for most other N-dimensional state spaces 
such a reduced representation is not feasible. In order to distin- 
guish between oscillators with different natural frequencies, we 
color the dots according to the standard color spectrum: the neu- 
rons correspond to the low end of the spectrum (close to -Jt)(red), 
neurons at the high end (close to jt) (blue), while those in between 
occupy the middle part of the spectrum (orange/yellow/green). 
To show how the system evolves from one instant to the next, 
we plot a series of snapshots of the system at different times 
(see Figures 11B-D, for example). This allows us to observe the 
behavior of individual neurons at the same time as we witness the 
collective evolution of the system toward an attractive state. 

3. RESULTS 

3.1. SINGLE NEURON BURST DYNAMICS 

We first examine the behavior of single neuron model 
Equation (2.1) as the applied input current I is brought close 
to the threshold for generating spikes or bursts. For the given 
parameters In Equation (2.1) a neuron is excitable. Figure 1 
depicts the relationship between applied input and parabolic 
bursting pattern. We are only interested in the behavior of this 
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FIGURE 1 | Shown here trace of membrane potential and calcium (A,C). In (B,D) membrane potential dynamics is shown for two 
dynamics. Fast spikes rides on a slow modulation of calcium. Slow different cases (A) / < l c and (B) I > l c in the single neuron 
subsystem moves Ca back and forth across SNIC bifurcation points model. 



system for low current values where the resting state of mem- 
brane voltage is sufficiently depolarized below —40 mV. When 
the applied input current / is below a critical value membrane 
potential V(t) maintains their steady state value and for val- 
ues greater than the threshold exhibits bursting behavior. For 
the parameterization used here we find that at I > 60 steady 
state destabilizes exhibiting multispikes. When the applied input 
current is further increased a neuron make transition from burst- 
ing to spiking behavior. In order to observe a typical spiking 
behavior we set I > 70. To get an intuitive understanding about 
the relationship between slow and fast subsystems, Rinzel and 
Lee analyzed this model by varying ca (a variable in the slow 
subsystem) as a bifurcation parameter to report that parabolic 
bursting is obtained from an oscillation in the slow subsystem 
that periodically moves the ca variable back and forth across the 
SNIC bifurcation, to link the steady state solution of this sys- 
tem to (quiescence state in Figure 1A) the branch of periodic 
solutions (Figure IB) and vice versa. Time series of fast variable 
shows that the interspike interval is relatively longer at the begin- 
ning and end of each burst. As has been shown by numerous 
authors oscillation for the fast dynamics is obtained when the 
slow variables are held fixed; it is where the saddle-node-loop 



bifurcation occurs. There is a clear threshold below which there 
is a unique stable fixed point. Parabolic bursting can occur with- 
out having any bistability in the spike generating process. One 
way to achieve parabolic bursting behavior without requiring any 
bistability in the generating process and moreover, mathemati- 
cally tractable would require a generic description like the one 
shown in Equation (2.2) (see section 2). From numerical results 
we find that as the applied input current I — > 2, time period 
T — > 00. Applied input current can be tuned such that it is pos- 
sible to obtain parabolic bursts of desired interburst gap. The 
time evolution of a single neuron activity is shown in Figure 2, 
where a membrane potential like variable V(f) = — cos[9(f)] 
is plotted by numerically integrating Equation (2.2). Temporal 
dynamics shows regular parabolic bursting behavior. For I < I c 
a neuron fires few spikes before it settles into a steady state. 
For I > I c (I = 2.01, n = 5) neurons exhibits parabolic burst- 
ing behavior. Based on the qualitative similarity in the burst 
pattern with parabolic bursting neurons (At the start and the 
end of the active phase the spike frequency is smaller com- 
pared to the middle of the active phase as can be seen in 
Figure 2 detailed model is substituted to investigate the network 
effects. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



March 2013 | Volume 7 | Article 20 | 6 



Roy and Jirsa 



Network time scale and dynamics 




500 1000 1500 2000 2500 3000 3500 

Time[ms] 




500 1000 1500 2000 2500 3000 3500 

Time [ms] 



FIGURE 2 | The temporal dynamics of the phase model of spike-burst neuron. V(t) = — cos[0J is plotted as a function of time, l c = 2.01 , n = 5, for 
(A) / < /„ and for (B) / > /„. 
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FIGURE 3 | Shown here trace of synaptic variable and approximated 
synaptic variable. Traces are plotted for a spiking regime of our network at 
/ = 80, this external current is applied to each neurons in this population. 
Panels (A-D) are generated for low to high fi synaptic time constant values, 
figure shows approximation breaks down progressively as we go to higher 
p values or access slower time scale. Approximation holds for faster time 
scale of oscillations. 



3.2. NETWORK DYNAMICS OF PARABOLIC BURSTING NEURONS WITH 
HETEROGENEITY 

To understand the influence of heterogeneity in the applied input 
current and the coupling strength in a network of single neurons 
exhibiting parabolic bursting we use Equation (5) and parameters 
as described in section 2.4. We use fast excitatory synapses to cou- 
ple these units. When the synaptic coupling is sufficiently fast, the 
coupling tends to push the neurons toward anti-synchrony (Wang 
and Rinzel, 1992; Friesen, 1994; Van Vreeswijk et al., 1994). 
Moreover, several studies have observed emergence of multistable 
solutions in their mean field network with parameter heterogene- 
ity (Assisi et al., 2005; Jirsa and Stefanescu, 20 1 1 ) . Our motivation 
is to go toward this particular direction to capture the rele- 
vant network dynamics at the population level. In particular to 
understand the combined effect of heterogeneity in the firing 
rate threshold (biophysical model) with the fast time scale of 
activation-deactivation of synapses in the coupling; the interplay 
between these two critical factors in spike/burst timing at the 
population level is largely unknown. In our formalism their indi- 
vidual and combined influence on the network dynamics become 
clearly visible. Typical time course of such responses of synap- 
tic variable in our model simulation are shown in Figure 2). Fast 
synapse approximation holds as long as the variable s,- relaxes 
much more rapidly than Vj, in which case we may apply a qua- 
sistatic approximation to reduce s; further in Equation (5), s, ~ 0 
allowing us to adiabatically eliminate (3, and set the synaptic 
variable via an approximation as s; = E . The mean 

l + p + exp^j 

synaptic action can be formulated as T = jj Ylf= l s i> where 
a s (Vj) = (i + exp (_y./2)) - The synaptic constant K is the same 
for all the neurons. Figures 3A-D shows kinetics of excitatory 
synaptic variable s, (plotted with black solid lines) for different 
f5 values. Mean synaptic variable (plotted with dotted lines) for 
the same set of values of time constant f5 shows dissimilar tempo- 
ral response compared to s; for higher time constant values. For 
smaller time constant values simulation provides relatively better 
aggrement as can be seen from Figure 3. We numerically integrate 
the above network to investigate how the mean population burst 
changes with time as a function of spread of applied stimulus A7 
and mean field coupling strength K. Firing patterns in this net- 
work are shown with array diagrams in Figures 4A,B. For small 
spread in the applied stimulus and sufficient coupling strength 
A7 = 0.001, K = 0.7 nearly burst synchronization takes place. 
Moreover, in the array diagram we detect clusters of synchronous 



states which fires in a wave-like pattern. Corresponding time 
series of mean quantities such as the membrane potential 
Vmean(f) i n Equation (11), mean field T(t) shows periodic activ- 
ity in Figure 4C. Membrane potential spiking activity is nearly 
synchronized across population of neurons in Figure 4. On the 
other hand, for the IN state mean membrane potential fast decays 
to zero and shows subthreshold fluctuations about mean zero. 
Response of mean membrane potential is more suppressed com- 
pared to their mean field oscillations between [0,1]. Amplitude 
of mean field T(f) changes in time systematically but fluctuates 
about the mean value of 0.5 instead of approaching zero val- 
ues as can be seen in Figure 4D. Population burst synchrony is 
observed for many different parameterization, for one such choice 
of parameter Al = 0.002, K = 0.8, an array diagram is computed 
and plotted in Figure 5A. As can be seen in the figure a wave- 
like spread of activity. In Figure 5B various time series plots of 
population burst synchrony is shown across 10 neurons. In order 
to identify different network states for all possible combination 
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FIGURE 4 | Array diagrams are shown in (A,B) for two distinct network 
states. In a nearly coherent states with clusters of synchronous bursting 
activity AI = 0.001, K = 0.7, in (B) incoherent states for A/= 0.12, K = 0.01. 
Nearly coherent states showing dynamical clustering effects and wave-like 
activity spread. Membrane potential time series is shown for all the neurons 
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exhibiting spiking dynamics both in the coherent and incoherent states. 
Mean membrane potential amplitude decreases and converge to a stationary 
solution. Standard deviation shown in (CD) shows growth in time. Mean 
field traces shows periodic activation and deactivation in the coherent state. 
In the incoherent state mean field amplitude systematically decrease in time. 



of two parameters K, AI we carry out a grid search and com- 
pute the values of R in Equation (12). Global order parameter 
measure identifies three distinct network states in the parameter 
space as shown in Figure 6A. For low coupling values K, order 
parameter shows fluctuations about mean zero. In this regime 
each neurons activity is mainly driven by their firing rate thresh- 
old and displays largely incoherence. For medium values of both 
coupling strength K and stimulus spread A7 network exhibits 
a hybrid state (some neurons are firing and some of them are 
silent). For very small values of stimulus spread and medium to 
high K values nearly burst synchrony appears. Temporal dynam- 
ics of membrane potential activity V(t) for four neurons are 
plotted in Figures 6B-D for three arbitrary parameterization of 
our network model. In Figures 6B,D PO state is shown where 



one neuron is spiking or bursting and three neurons are silent. In 
Figure 6C all neurons are showing nearly synchronized parabolic 
bursting behavior. 

3.3. NETWORK DYNAMICS OF PARABOLIC BURSTING PHASE MODEL 
WITH HETEROGENEITY 

In this section, we use a phase network with mean field cou- 
pling to get some insights about the novel network states observed 
in (K, AI) the parameter space of the full network model. 
Coupling between each phase neuron via a mean field is formu- 
lated in section 2.4. Numerically we integrate Equation (6) to 
compute time averaged membrane potential, mean field T as in 
Equation (10) (see section 2), global measure of coherence Rq as 
a function of K, AI a parameter combination which is used in 
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FIGURE 5 | In (A) array diagram showing firing pattern in a population of 
100 neurons. Only 10 neuron index are shown for clarity. Horizontal axis 
is always time and vertical axis is labeled as neuron index. Green color 
corresponds to no firing activity or quiescence. Purple pixels corresponds to 



parabolic bursting activity of each individual neurons which are locked in time. 
In (B) time series data for membrane potential of V(t), Vmean(f), T(f), and std 
Knean(f) are plotted for 10 neurons. Mean population burst synchronizes in 



the detailed network model. Time evolution of the above quan- 
tiles are shown in Figure 7 for a parameterization K = 0.8, AI = 
0.001. The parameter choice is the same as the full network inves- 
tigations. With this combination of parameters all the neurons 
synchronously spikes. Mean membrane potential-like quantity 
Vmean(f) oscillates in phase with synchronized spike activity as 
plotted in Figure 7B. Here, n a quantity which determines the 
number of spikes per burst is kept at n = 1. Mean field F also 
shows up and down states (Locked in time) and act as an oscil- 
lating drive to each individual neurons. The time series of the 
global order parameter -Re(f) for synchronized spiking is periodic 
in Figure 7E. Next, we show in Figure 8 temporal evolution of the 
mean quantities for the choice ofK = 0.8, AI = 0.5. For medium 
values of mean field coupling strength and stimulus spread net- 
work shows PO behavior, where some of the neurons are firing 
incoherently and others are completely silent. This means for 
some parameterization network has two stable branches of solu- 
tions. It is important to note PO state of the network was 
observed in the full network for a comparable parameterization 
(see Figure 6). Time series for 10 neurons and their order param- 
eter evolution in time is plotted in Figure 8. Three neurons are 
completely silent while other seven neurons are bursting with 
variable inter-burst intervals. As there is no noise in this system 



and coupling magnitude is set at high values as in the case of sync, 
this variability must be introduced by the heterogeneity in their 
individual firing rate threshold via stimulus spread. 

Figure 6 shows the parameter space diagram for the full and 
phase models presented in Equations (5) and (6-9). Phase bound- 
aries are calculated by computing the mean field for both full 
and the phase model for different combination of (K, I) values 
on a two-dimensional grid. In the following subsection we would 
lay out the details for obtaining the phase transition boundaries 
semi analytically. Interestingly, Over a wide range of (K — AI) 
values the collective dynamics of the two networks primarily 
show three distinct regions of interest which are close to each 
other in the parameter space. For sufficiently large K values hold- 
ing A7 fixed to a narrow range of values near zero, the system 
converges to a state of partial oscillations in which the some of 
the neurons are not firing altogether, while the others display 
IN oscillations. Very large K values result in damping of oscil- 
lation activity and all the neurons stops firing altogether. The 
stability state of locking is much more difficult to achieve and in- 
fact we found distinct branches in their rotation number, these 
states should all be regarded as variants of 1:1 locking, and there- 
fore we lump them together in the locked region of the stability 
diagram. With further increase in AI, parameter heterogeneity, 
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FIGURE 6 | Phase diagram of mean synaptic action variable is shown 
as a function of 2D parameter space of stimulus spread A/ and 
excitatory coupling strength K. In the partial burst regime labeled as 
PO in (A), a subset of neurons are not firing at all as their respective 
drives are below their individual firing thresholds. Heterogeneous 
dynamics between synchronized population spiking activity and oscillation 
frequency death response for PO state is displayed in (B). Nearly 



0 2000 4000 6000 8000 10000 
Time [ms] 

synchronized population of spike/burst activity lumped in a regime labeled 
as SR [corresponding time series is displayed in (C)] and incoherent 
population spike/burst activity is lumped into a regime called IN. In the 
incoherent regime mean field values stay close to zero with substantial 
subthreshold fluctuations. In (D) multi stability of PO state is displayed 
again; now between population burst and fixed point dynamics for an 
entirely different parameterization. 



successively more neurons peel away until eventually the entire 
population is IN. 

4. PHASE DIAGRAM USING SEMI ANALYTICAL METHODS 
FOR MEAN FIELD PHASE MODEL 

Mean field coupled neurons in phase model is described in 
Equation (6). Let's rewrite the mean field equation explicitly. 



9; = (F(9,) - T sin9,(cos9; - Vth)) 



(15) 



where r - K ^ N 



EN 
j = 



^l+(S + exp^-- 

In a semianalytical approach we would like to understand the 
phase transitions between three distinct network states discov- 
ered in two networks. For the IN states where the average firing 
frequency increases monotonically plotted in Figure 10, the 6; 
are all distributed across the closed orbit in a unit circle. This 
leads to the following phase evolution equation SR state may 
undergo instability either through parameter changes of K or 
AJ and make phase transition to either IN or PO state. Mean 
field T approaches a stationary density as the number of neurons 
are increased in both PO and IN state (see Figure 4D). Hence, 



T approaches some positive real number for these two states. 
When varying K, we consider small perturbations |x to the 
SR solution 6 = 6, = 0. With 6 = 9, ■ = 0 + (jl Equation (15) 
becomes 9, = (i = F(|x) + j sin (2[i) and linearization yields 
9; = (f'(0) + V)[L. Moreover, SR state may gets phase locked at 
9 = 7T (subpopulation clusters). Hence, 9 = it may get destabi- 
lized as we changed the width of heterogeneity by changing AJ 
or the coupling strength K. Similarly, we consider small per- 
turbations |jl about solution 9 = n. Hence, we can write 9; = 
Jt + |x, 9,- = (1 = _F(tt + jjl) + j sin(27t + 2[i) and linearization 
yields 9, = (F'(n) + F)\i. With F(9) « I - cos 9 - cos £ for the 
SR state, we find that F'(Q) ~ — sin(9) — - sin - will be gener- 
ally small for 9 = 0, Tt. SR state solution hence becomes unstable 
when F'(0, tt) + T = 0, which suggest almost a vertical critical 
line between SR and PO, IN state. The bifurcation route from PO 
(multistable state) to IN solutions as the parameter AI increases 
is less conclusive in the framework of the circular approximation, 
since in the previous stability analysis the only /-dependent term 
is F'(0, Tt), which is very small, hence higher orders of the approx- 
imation must be considered. We use the following ansatz: If r is 
the radius of a unit circle, any smooth deformation from a unit 
circle can be approximated as, r(6,) = 1 + €/j(9,). Hence we can 
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FIGURE 7 | Time series of (A) V{t) = - cos9 ; (t), (B) V mea „{t), (C) std V mea „{t), (D) mean field r, and (E) order parameter R„{t) are plotted for 10 
neurons. All neurons are spiking in synchrony and time locked. The parameter values are K = 0.8, Al = 0.001 . 



compute the non-linear flow contribution with the above first 
order correction term as F(Q,) + eH(0j). It is possible to explic- 
itly determine 77(9,) for a certain choice of h(Q,) and moreover, 
77(9;) has a periodicity of n, that is H(9 + it) = H(0). Thus the 
linear stability analysis about the fixed point 9* = 0 + \i gives 

ji = (F'(0) + € H'(0) + r)(i = (F'(jt) + €H'(it) + r)|x (16) 

From the above equation with F'(0) ~ 0 and the jt-periodicity 
of 77(9,), we find that the two fixed points at 0, tt lose stability 
at the same time for increasing A7 and as a result leads directly 
to the IN state. Since 77' (9,) ~ 7, scales linearly for fixed u,, we 
can also estimate the critical line of transition in the parame- 
ter space in Figure 11 which separates PO state from IN state. 
For the critical line: 77'(9,) = m(I — I C )P where m is the slope 
of this line and m > 0 allows for destabilization. Hence the 
critical condition is €77' (0) + T = 0. By substituting the depen- 
dence of 77' (9;) on (7, 7 C ) and in turn dependence on A7 one 
can write em(7 — I c ) p + F = 0. This implies coupling strength 
K = — «(A7 + AI C )P for (m > 0) andp is some exponent repre- 
senting a scaling relationship near saddle-node bifurcation. Thus 
the critical condition is \K\ = em(A7 — A7 C ), which serves as a 
convenient guide to numerically compute the stability line sep- 
arating PO region from IN. Next we try to obtain analytically 



the stability boundary between INC and PO oscillation states in 
the infinite-N limit, it turns out that the IN and partial oscilla- 
tion states can be made steady in our system. The possibility of 
doing so was suggested by the numerical results. In numerics we 
observed that as the number of neurons N is increased, the order 
parameter R t approaches a constant for both these states Figure 8 
and the oscillators tend to arrange themselves in a stationary 
distribution around the circle Figure 11. The way to approach 
these two states analytically, therefore, is to first write down the 
appropriate infinite-N analog of our model. 



— f — 
df + 39 



!/ / pin pl+ AI 
m ~\\L Si ai r P@,t,r) g (i') l u'dff 



x sin(29) 



0 



(17) 



The above equation is the infinite-N analog of continuity 
equation introduced earlier. It is a non-linear partial integro- 
differential equation for the number density /(9, f, go). In addi- 
tion we demand / to be non-negative, 2tt periodic in 9, and we 
impose the normalization 



(•2it 

Jo ] 



f(Q, t,I)dQ=l, 



(18) 
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FIGURE 8 | Temporal evolution V/(t) for N = 10 neurons are shown here 
for an arbitrary parameterization K = 0.8, A/ = 0.5. K value is unchanged 
from previous figure. Stimulus spread A/ is changed. Time series of order 

For incoherence and partial oscillation the above system 
tends toward a stationary distribution of phases in time. 
The above two states are the fixed points of the station- 
ary density in the continuum limit. To solve for the fixed 
points we set 4/ = 0 in Equation (10). let's assume that 
/o(6, w) be the stationary phase density and Vo = [F (0) — 
((fo* fi-Ai r f(P> t,I')g(I')dI'dQ')) sin(2Q)] be the velocity 
field. Then one can write 



parameter Re(t) undergoes statistical fluctuations of magnitude 0( 
about some positive constant value. After 8000 time points, dynamics is 
truncated assuming network dynamics settles into a steady state. 



where L(7) is a constant which is determined exactly by 
using normalization condition. Depending on it's applied drive 
I, neuron's steady state behavior falls in the following two 
categories: 



Case (i) When I < < F implies 



90 



(fovo) = 0 =► / 0 v 0 = 1(7) 



(19) 



vo(9, 1) = F(9) - T sin(26) = 0 



(20) 
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FIGURE 9 | In (A,B) an array diagram is shown for phase network model 
for a parameter combination. In (A) K = 0.61 , AI = 0.001 , n = 1 spikes 
only and in (B) K = 0.001 , A/ = 0.15, n = 3 bursts only. Almost always, near 
synchronous burst states are observed for high K and low AI values. In (C) 
corresponding time evolution of Vj(t) is shown for all 10 bursting neurons. In 
(D) temporal response of Vmean is shown for the coherent state of our 



network. In (E) standard deviations of V mem is plotted as a function of time. 
In (F) mean field r vs. time for the coherent state is shown. (G) displays 
temporal dynamics of order parameter. Average firing frequency as described 
in Equation (14) is plotted in (H) for the parameter combination of K = 0.001, 
A/ = 0.15. Panel (H) further demonstrates phase locking behavior among all 
the neurons. 



Case (ii) When I > > F neuron fires incoherently and typically 
individual phases follows an uniform distribution about the unit 
circle. In this case the velocity field turns out to be, 



v o (0,I) = F(G) - rsin(26) 



(21) 



Fixed point solution demands that the density must be inversely 
proportional to the velocity: 



/o(e,J) 



F(6) - T sin(29) 



(22) 



In the IN state, neurons driven by different external drives are fir- 
ing at different phases, however, their collective state is close to 
being stationary. Every neuron belong to Case (ii) as described 
above. Further, it is possible to derive nearly an exact relation- 
ship between K, AI that gives the transition from case (i) to 
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case (ii) as described above. As shown before, in case of a finite 
size network such a relationship in the first order perturbation 
em(7 — I C )P + T = 0 does exist. In this scenario those neurons 
with a minimum bound on their applied drive 7 m i n reach cessa- 
tion of firing as we find from numerical simulations. They then 
fall into the Case (i) above where mean field T exerts much bigger 
influence on the dynamics and overall effect is damping of fir- 
ing activity. The first neurons to stop firing are the ones which do 
not cross the threshold for firing which in this case I > 2. Then 
the boundary that separates IN from PO in the phase diagram is 
almost a straight line given by, 

\K\ = €ot(A/ - AI C ) (23) 



Hence, both finite and infinite analog of our network iden- 
tifies the putative transition boundary between IN and PO 
states. Now from numerical simulations we find Andronov- 
Hopf (AH) bifurcations leads to the transition from INC to 
SR solutions in the Figure 6 near K, I values close to zero. 
It is equivalent to look at the imaginary eigensolutions that 
arise due to the instability of the IN state. This instability 
requires calculation of higher order perturbation terms of the 
stationary density obtained at the IN state of our network. 
This is out of the scope of our paper, however, we show a 
numerical fitting result which gives an empirical relationship 
between K and A7 to quantify the transition boundary between 
IN and SR states. Assuming 6 is the perturbation to the IN 
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frequency of individual neurons are shown in a unit circle. Each position in a 
circle corresponds to a particular phase and color coded according to the 
scheme described in section 2.8. 



solution we can express a relationship between K and Al as 
follows, 

\K\ = a 0 e + a^ 2 + a 2 e 3 + 0(« 4 ) (24) 

Equation (24) gives us an empirical relationship between param- 
eters upto fourth order perturbations for the bifurcation of a 
limit cycle. Optimization of the above equation gives coefficients 
a 0 = I", fli = 0, fl2 = respectively. Next, we substitute the 
amount of dispersion Al into the perturbative term € to obtain 
the boundary between IN and SR state. Taken together we can 
write, 

128 



1*1 



-AI + —-AP + 0(AP) 



(25) 



Results are shown in Figure 12 in the (K, Al) plane using 
Equations (25) and 23. Critical lines obtained semi analytically 
qualitatively agrees well with the numerical results that cap- 
tures various network states in both these models with purely 



excitatory coupling. In Appendix, we show a stability calcu- 
lation for an inhibitory coupled mean field network in the 
infinite analog limit. From numerical simulations we find that 
the results are independent of the number of spikes n per 
burst. 

5. DISCUSSION 

One of the most frequent assumption employed in simulations of 
large neural networks is that the whole network can be lumped 
into small aggregates of collective unit (sometimes called a "neu- 
rocomputational unit") exhibit a sufficiently similar dynamical 
behavior. Consequently, the network that instantiates this ensem- 
ble, consisting of thousands of excitatory and inhibitory neurons, 
it is considered to display a synchronized behavior with no other 
significant temporal features for the dynamics of the large scale 
network. The main reason for this assumption, is the imprac- 
tical large computational time arising from too many details 
considered in the large network properties. In this paper, we 
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have analyzed the behavior of a neural network that serves as 
a good example of such a unit, namely a mean field coupled 
bursting ensemble. First, we have investigated a Hodgkin-Huxley 
type detailed biophysical model widely employed in theoreti- 
cal and computational neuroscience with global coupling. We 
found that the dynamical features of the network are far more 
complex then the ones corresponding to synchronized or rest 
state behavior. The network dynamics depends critical on the 
balance between firing rate threshold dispersion and mean field 
synaptic coupling strength; in fact, the synchronized state can be 
found only for a specific range of parameters typically involving a 
large or medium values for the coupling strength and low val- 
ues of dispersion. On the other hand, for large dispersion and 
weak coupling strength values both networks display purely IN 
behavior. In the IN state, individual neurons are driven by differ- 
ent external drives results in firing at different phases, however, 
their collective state is close to being stationary. This stationarity 
in the density distribution led us to formulate scaling relation- 
ship between coupling strength and dispersion parameter. One 
interesting finding is that, when mean field exerts a greater influ- 
ence than parameter dispersion; it causes shutting down of the 
neural activity in some neurons. In this parameter range, we 
find interesting dynamical behavior such as partial activity. In 
order to address the problem of the high computational cost 
of such an implementation, we have further developed a self- 
consistent mathematically tractable mean field coupled phase 
model following (Assisi et al., 2005; Ghosh et al., 2009; Jirsa and 
Stefanescu, 2011), but incorporating a higher degree of realism. 
Rather than finding the most appropriate type and number of 
dimensions that could minimize certain error functions or cap- 
ture statistical variance in the full network, we have focused our 
attention on understanding a phenomenological burst genera- 
tion model system which captures the most important network 
dynamics of bursting units at the population level. Collective 
activity of synapses is described by a mean field which relies on 
instantaneous rise and decay time (Roy et al., 2011). This mean 
field is then employed in the coupling to individual neurons 



to describe phase network. Together, we investigate this popu- 
lation of neurons coupled to a common mean field drive and 
heterogeneity in their threshold for spikes/bursts. Our detailed 
analysis demonstrated that the reduced representation manages 
to recreate correctly the topology of the mean field amplitudes 
of the original system for various parameter scenarios. In the 
full network, In the thermodynamic limit (N oo), a collective 
state becomes coherent if 5V me an(0 = [V me an(0 - V mean (t)] is 
non-stationary (i.e., an oscillating global potential V me an appears 
for a coherent case) and also, the correlated mean field F(t) 
activity appears oscillatory. In the phase network, global order 
parameter is computed by averaging the contributions of all 
microscopic spikes within a burst in order to obtain a simi- 
lar degree of ordering of spikes/bursts as in the full model for 
identical parameterization. Hence, for a dynamical behavior such 
as coherence-incoherence transition macroscopic order param- 
eter gives us a crude approximation of burst timing. From a 
more general perspective, despite its limitations because of the 
consideration of purely excitatory or inhibitory network, it can 
be readily extended to study networks with mixed coupling. 
Moreover, the analytical approach to estimate the scaling rela- 
tionship and transition boundaries between the IN-PO-SR states 
is not restricted to small scale network only. With global cou- 
pling, as the size of the network grows the boundaries may shift 
leading to a different parameterization than the one displayed 
here; however, underlying bifurcations remain the same. We have 
demonstrated this in our work by analytically deriving a low 
dimensional mean field amplitude reduction for a inhibitory cou- 
pled mean field network in the continuum limit. In this case, 
all the relevant dynamics of an infinite dimensional network 
in Equations (29) and (30) is captured by a two dimensional 
representation of the reduced mean field population given by 
Equation (40). Thus, using this approach, we derive analytically 
a low dimensional representation of the network dynamics and 
we show that the main features of the neural population's col- 
lective behavior can be captured well by the dynamics of a few 
cortical nodes exhibiting spiking as well as bursting behavior. 
While it is true that strong reductionist assumptions are common 
(sacrificing many of the biological realism of a network node's 
dynamics) in large-scale network modeling, these assumptions 
are usually made ad-hoc on the network node's dynamics and 
limit the network dynamics to a small range. We emphasize here 
that because of the "near to synchrony" assumption, neural mass 
models cannot capture complex dynamical features such as multi- 
clustering, oscillator death or multi-time scale synchronization. 
Evidently a reduced small scale network model is desirable to 
serve as a node in a large scale network simulation whereby dis- 
playing a sufficiently rich dynamic repertoire. Here it is of less 
importance to find a quantitatively precise reduced description 
of a neural population; rather more importantly, we seek a com- 
putationally inexpensive population model (this means typically 
low-dimensional) which is able to display the major qualitative 
dynamic behaviors (synchronization, rest state, multi-clustering, 
etc.) for realistic parameter ranges as observed in the total pop- 
ulation of neurons. Our approach may offer a viable alternative 
to the neural mass models currently used in the literature. By 
comparison, our model offers the possibility to account for such 
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features (temporal details of their spiking activity considered 
irrelevant for the dynamics of the large network) at a very low 
computational cost. Therefore, the type of reduced representa- 
tion discussed in this paper qualifies as a good candidate for a 
"neural unit" in computational simulations of large scale neural 
networks. 
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APPENDIX 

MEAN FIELD REDUCTION FOR INHIBITORY SYNAPTIC COUPLING 

Here, we extend our network in the continuum limit in the 
presence of inhibitory coupling. V t h is held negative. (N — > oo), 
where the state of the coupled system can be described by a den- 
sity function /(9, 1, t), where / is defined such that the fraction 
of neurons with phases lying between 6 and dQ and applied drive 
between I and dl is given by/(6, J, t)dQdI (Antonsen et al, 2008; 
Ott and Antonsen, 2009). The applied stimulus are drawn from a 
distribution g(I) such that 



/oo pin 
/ f(Q, I, t)dQdI = 1 
-oo Jo 



f*f{B,I, t)dQ = g(I) 
Jo 



(26) 
(27) 



Using the above it is easy to see that the expression for velocity 
becomes 



v(Q,I, f)=I+j. + |(z2 + z1*)+f)< 
- (l + ^(zl+z2*)+l)e' e ] 



jO/n p —ito/n 



2i 2i 



(33) 



* indicates the complex conjugate. The distribution function can 
be expressed as a Fourier series 



mi, t) 



2 it 



!+£>(!, t)e M + c.c. 



k=l 



(34) 



For the conservation of currents I the continuity equation is 
written as 



9/ 9(/v) =Q 
dt 39 



(28) 



The above infinite dimensional system is difficult to analyze. 
However, the "amazing" anstaz of Ott and Antonsen (2009) has 
been shown to be successful in obtaining the low-dimensional 
description of the globally coupled phase oscillators. The anstaz 
impose a restriction on the fourier coefficients: 



In order to make the coupling amenable to analytical study we 
use a pulse-like function for the mean field r = fll + fol(l + 
cosO). Response to the mean field by individual neuron's R(Q) = 
V t hs'm(Q), containing only single Fourier component, a choice 
motivated primarily due to the tractability of the resulting model. 
Further, Vth = — 1 f° r the convenience of calculations without 
losing any generality of our results. The velocity v(9, 7, t) in 
Equation (28) is now written as 



/oo t>2it 
/ (l + cos(9))(- sin 9) /(§,/, t)dM 
-oo Jo 



-oo J0 

Fsin(9)-Fsin(9/n) 



(29) 



f k (i t) = oki, or 



(35) 



for k > 1 and has been shown to be a reasonable guess under 
different scenariors (Ott and Antonsen, 2009). This restricted 
class of functions readily reduces our continuity equation to an 
9-independent form 



dt 



2 V 2 / 



i0 



1 + 1/k 



^1-1/" 



(36) 



where without loss of any generality we are using sin functions 
instead of cos functions in Equation (6). 

In the continuum limit the order parameter z can be 
defined as 



with zl satisfying 



/oo 
Ajr*(I, t)g(l)dl. 
-OO 



(37) 



z(t) 



n 



2 it ,J0 



(e m + e- M ) 



f(Q, I, t)dQdI 



(30) 



It's a linear sum of two complex order parameters and one could 
in principle unfold the entire dynamics of the network in any one 
of the manifold given above. Here, 



/oo /*2jt 
/ e B f(Q, I, f)dUl 
-oo J0 

/oo />2jt 
/ 
-oo Jo 



e~ ,9 /(6, l t)dQdI 



(31) 
(32) 



If we assume that g(7) is a Lorentzian distribution function 



Tt[(7 - 7 0 ) 2 + 1] 



(38) 



z(t) can be evaluated by contour integration with poles at 7 = 
7o — i and we obtain the exact evolution equation of order 
parameter z 



dzl 
dt 



HqzI — zl + 
< gl l + l/n 



1 + fzl + F l+|zl*+F 



zl 



l - l/n 



(39) 
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The above equation can be expressed in polar coordinates if 
we substitute zl = p\exp(i<\>\) giving evolution equations for 
Pi and 4>i 



dt 



dt 



rPi(l-Pi) 



F , 



2 P \y 



+ 



Fp 



1 / -l/n ,, , s 

-\Pi - Pi )cos(4>i/n) 



1 

>i + - COS((J)i) 



(40) 



= k 



Pi + 



Pi 



sin 4>i 



Pi 



sin(<|)i) 



l/n . 
Pi + 



Pi 



l/n 



sin((t>i/n). 



(41) 



For the Lorentzian distribution function the above equation is 
exact, However, we do not find any deviation of the above results 
for any other unimodal distributions of our firing threshold 
(such as uniform distribution) such as the one considered in the 
numerical simulations with excitatory coupling. The above two 
dimensional system can be solved numerically to identify the full 
network states and the corresponding transition boundaries. 
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